% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
% Generate firm panel

I               = 2000; % number of firms to start with
T_sim           = 5000; % number of periods
delx            = 300; % periods to delete

rng('default');
rng(1);
panel_b_grid        = zeros(I,T_sim); % debt (grid point)
panel_z_rank        = zeros(I,T_sim); % grid point of Z
stat_z              = sum(dist)/sum(dist(:)); % distribution over productivity
panel_z_rank(:,1)   = randsample(1:lenz,I,true,stat_z); % initial Z vector (grid)
for xx = 1:I    
    panel_b_grid(xx,1)   = randsample(1:lenb2,1,true,dist(:,panel_z_rank(xx,1))); % initial B vector (grid)
end

% Policy functions:
panel_b             = zeros(I,T_sim); % LT debt choice
panel_b_iss         = zeros(I,T_sim); % LT debt issuance
panel_bs            = zeros(I,T_sim); % ST debt choice / issuance
panel_q             = zeros(I,T_sim); % assets
panel_k             = zeros(I,T_sim); % choice of capital
panel_e             = zeros(I,T_sim); % choice of etilde
panel_default       = zeros(I,T_sim); % default decision (yes/no)
panel_epsbar        = zeros(I,T_sim); % choice of default threshold
panel_iwspread      = zeros(I,T_sim); % issuance-weighted spread
panel_ltspread      = zeros(I,T_sim); % LTD spread
panel_exit          = zeros(I,T_sim); % indicator: firm exits end of period
panel_alive         = ones(I,T_sim); % indicator: firm is alive
panel_age           = ones(I,T_sim); % firm age
panel_ic            = zeros(I,T_sim); % debt issuance costs
panel_y             = zeros(I,T_sim); % output
iwspread_pol        = reshape(ltspr,size(ltd_pol2));
ltdspread_pol       = reshape(ltspr,size(ltd_pol2));

for t = 1:T_sim
    idx                 = (panel_z_rank(:,t)-1)*lenb2 + panel_b_grid(:,t);
    panel_b(:,t)        = ltd_pol2(idx); % choice of LTD
    panel_b_iss(:,t)    = panel_b(:,t) - debtgrids(idx); % LTD issuance
    futuredebt_stay     = bp_pol_stay(idx);
    futuredebt_up       = bp_pol_up(idx);
    panel_bs(:,t)       = std_pol2(idx); % choice of LTD
    panel_k(:,t)        = cap_pol2(idx); % choice of capital
    panel_e(:,t)        = eq_pol2(idx); % choice of etilde
    panel_epsbar(:,t)   = epsbar_pol2(idx); % choice of default threshold
    panel_iwspread(:,t) = iwspread_pol(idx); % spread
    panel_ltspread(:,t) = ltdspread_pol(idx); % ltd spread
    panel_y(:,t)        = PARMS.Z(panel_z_rank(:,t)).*(panel_k(:,t).^PARMS.PSI .* lab_pol2(idx).^(1-PARMS.PSI)).^PARMS.ZETA; % output
    panel_ic(:,t)       = PARMS.ETA*(panel_bs(:,t) + max(panel_b(:,t) - debtgrids(idx),0)).^2;  % debt issuance costs
        
    % UPDATE THE PANEL: Draw new Z shocks.
    panel_z_rank(:,t+1) = panel_z_rank(:,t); 
    panel_b_grid(:,t+1) = futuredebt_stay;   
    n                   = randsample(1:I,round((1-PARMS.rhoz)*I));
    panel_z_rank(n,t+1) = min(lenz,panel_z_rank(n,t+1) + 1);
    panel_b_grid(n,t+1) = futuredebt_up(n); 
    epsshock            = 0 + PARMS.sigmaeps*randn(I,1);
    panel_default(:,t)  = epsshock < panel_epsbar(:,t); 
    panel_z_rank(panel_default(:,t)==1,t+1) = 1;
    panel_b_grid(panel_default(:,t)==1,t+1) = 1;
    
    part_nd             = (1-panel_default(:,t)) .* (panel_k(:,t) -  panel_bs(:,t) - PARMS.GAMMA*panel_b(:,t) + ...
        (1-PARMS.TAU)*( panel_y(:,t) + (epsshock - PARMS.DELTA).*panel_k(:,t) - PARMS.W*lab_pol2(idx) - PARMS.F - ...
        PARMS.C*(panel_bs(:,t) + panel_b(:,t)) ) );
    part_d              = 0; 
    panel_q(:,t)        = part_nd + part_d;
       
    exitshock           = rand(I,1);
    panel_exit(:,t)     = exitshock < PARMS.x;
    panel_z_rank(panel_exit(:,t)==1,t+1) = 1;
    panel_b_grid(panel_exit(:,t)==1,t+1) = 1;   
    panel_age(:,t+1)    = panel_age(:,t) + 1;   
    panel_age(panel_exit(:,t)==1,t+1)    = 1;
    panel_age(panel_default(:,t)==1,t+1) = 1;
    
end

panel_z_rank(:,T_sim+1)     = [];
panel_age(:,T_sim+1)        = [];
fprintf('Panel simulation done\n');

panel_b_grid(:,1:delx)  = [];
panel_z_rank(:,1:delx)  = [];
panel_b(:,1:delx)       = [];
panel_b_iss(:,1:delx)   = [];
panel_bs(:,1:delx)      = [];
panel_k(:,1:delx)       = [];
panel_e(:,1:delx)       = [];
panel_default(:,1:delx) = [];
panel_exit(:,1:delx)    = [];
panel_epsbar(:,1:delx)  = [];
panel_alive(:,1:delx)   = [];
panel_age(:,1:delx)     = [];
panel_q(:,1:delx)       = [];
panel_ic(:,1:delx)      = [];
panel_y(:,1:delx)       = [];
panel_iwspread(:,1:delx)= [];
panel_ltspread(:,1:delx)= [];
panel_z                 = PARMS.Z(panel_z_rank);

% Compute LTD share, leverage, and LTD share of debt issuance 
panel_maturity          = panel_b ./ (panel_b + panel_bs) * ((1-PARMS.GAMMA)/(1+PARMS.R))^4;
panel_leverage          = (panel_b + panel_bs) ./ panel_k;
panel_ltdshare_iss      = panel_b_iss ./ (panel_b_iss + panel_bs);

% Aggregate investment rate
K_agg                   = mean(panel_k);
I_agg                   = [0 (K_agg(2:end) - (1-PARMS.DELTA)*K_agg(1:end-1)) ./ (K_agg(1:end-1))];

panel_out               = min(1,panel_exit + panel_default);
idx_i                   = 0; % indicator for firms
P_aggK                  = cell(1,1); % Aggregate capital in each period of a firm's life
P_k                     = cell(1,1); % Panel for firm capital choice 
P_T                     = cell(1,1); % Panel for aggregate time period
P_std                   = cell(1,1); % Panel for firm STD choice (end of period)
P_ltd                   = cell(1,1); % Panel for firm LTD choice (end of period)
P_spread                = cell(1,1); % Panel for LT credit spread

panel_irate             = zeros(size(panel_k)); % Investment rate
panel_eqiss             = zeros(size(panel_k));
panel_diss              = zeros(size(panel_k));
panel_eqshare           = zeros(size(panel_k));
panel_xshare            = zeros(size(panel_k));
panel_i                 = zeros(size(panel_k)); % Investment (level)
panel_k_lag             = zeros(size(panel_k)); % Lagged capital
panel_y_lag             = zeros(size(panel_k)); % Lagged output
panel_deq               = zeros(size(panel_k)); % Equity financing
panel_dd                = zeros(size(panel_k)); % Debt financing

for i = 1:I % go through all firms   
    done_i = sum(panel_out(i,:));
    idx_t1 = 0; % reset time indicator to zero
    while done_i > 0
        % Update firm counter and time index
        idx_i   = idx_i + 1;
        idx_t   = idx_t1 + 1;
        idx_t1 = find(panel_out(i,idx_t:end),1,'first');       
        idx_t1 = idx_t - 1 + idx_t1;
        
        % Fill the panel for firm idx until it exits the panel
        P_k{idx_i}      = panel_k(i,idx_t:idx_t1);
        P_std{idx_i}    = panel_bs(i,idx_t:idx_t1);
        P_ltd{idx_i}    = panel_b(i,idx_t:idx_t1);
        P_T{idx_i}      = idx_t:idx_t1; % time
        P_spread{idx_i} = panel_ltspread(i,idx_t:idx_t1);
        P_aggK{idx_i}   = K_agg(idx_t:idx_t1);
        
        panel_irate(i,idx_t:idx_t1) = (P_k{idx_i}(1:end) - (1-PARMS.DELTA)*[0 P_k{idx_i}(1:end-1)]) ./ ([0 P_k{idx_i}(1:end-1)]);
        tmpq        = panel_q(i,idx_t:idx_t1); 
        tmpe        = panel_e(i,idx_t:idx_t1); 
        tmpeiss     = tmpe - [0 tmpq(1:end-1)]; 
        tmpdiss     = panel_b_iss(i,idx_t:idx_t1) + panel_bs(i,idx_t:idx_t1); % Debt issuance
        % Equity issuance divided by lagged capital:
        panel_eqiss(i,idx_t:idx_t1) = [NaN tmpeiss(2:end) ./ P_k{idx_i}(1:end-1)];
        % Debt issuance divided by lagged capital:
        panel_diss(i,idx_t:idx_t1)  = [NaN tmpdiss(2:end) ./ P_k{idx_i}(1:end-1)];
        % Equity share of external financing
        panel_eqshare(i,idx_t:idx_t1)= tmpeiss ./ abs(tmpeiss + tmpdiss);
        % External finance as fraction of total assets
        panel_xshare(i,idx_t:idx_t1)= [NaN abs(tmpeiss(2:end) + tmpdiss(2:end)) ./ P_k{idx_i}(1:end-1)];
        panel_i(i,idx_t:idx_t1)     = (P_k{idx_i}(1:end) - (1-PARMS.DELTA)*[0 P_k{idx_i}(1:end-1)]);       
        panel_deq(i,idx_t:idx_t1)   = tmpeiss;
        panel_dd(i,idx_t:idx_t1)    = tmpdiss;        
        panel_k_lag(i,idx_t:idx_t1) = [0 P_k{idx_i}(1:end-1)];        
        panel_y_lag(i,idx_t:idx_t1) = [0 panel_y(i,idx_t:idx_t1-1)];
                
        done_i = done_i - 1;
    end
    idx_i           = idx_i + 1;
    idx_t           = idx_t1 + 1;
    P_k{idx_i}      = panel_k(i,idx_t:end);
    P_std{idx_i}    = panel_bs(i,idx_t:end);
    P_ltd{idx_i}    = panel_b(i,idx_t:end);
    P_T{idx_i}      = idx_t:numel(panel_k(i,:));
    P_aggK{idx_i}   = K_agg(idx_t:end);
    P_spread{idx_i} = panel_ltspread(i,idx_t:end);

end
panel_irate(isinf(panel_irate)) = NaN;
numfirms    = size(P_k,2);
numpers      = cellfun('size',P_k,2);

%% COMPUTE MOMENTS FROM THE PANEL
ir_mean         = NaN(numfirms,1);
ir_mean_net     = NaN(numfirms,1);
ir_stds         = NaN(numfirms,1);
n               = 5;
for ii = 1:numfirms
    if numpers(ii) >= n      
        tmp             = (P_k{ii}(2:end) - (1-PARMS.DELTA)*P_k{ii}(1:end-1)) ./ (P_k{ii}(1:end-1));
        ir_mean(ii)     = mean(tmp);
        ir_stds(ii)     = std(tmp);
        tmp_agg         = (P_aggK{ii}(2:end) - (1-PARMS.DELTA)*P_aggK{ii}(1:end-1)) ./ (P_aggK{ii}(1:end-1));
        ir_mean_net(ii) = mean(tmp - tmp_agg);
    end
end

%% Export panel to compute correlations with size and age in Stata (Table 4)

panel       = zeros(1, 6); % id, time, k choice, std choice, ltd choice 
idx_t       = 0;
numfirms_s  = 5000;
for i = 1:numfirms_s    
    add                           = numel(P_T{i});
    panel(idx_t+1:idx_t+add,1)    = i; % id
    panel(idx_t+1:idx_t+add,2)    = P_T{i}; % time
    panel(idx_t+1:idx_t+add,3)    = P_k{i}; % capital choice
    panel(idx_t+1:idx_t+add,4)    = max(0,P_std{i}); % std choice
    panel(idx_t+1:idx_t+add,5)    = P_ltd{i}; % ltd choice
    panel(idx_t+1:idx_t+add,6)    = P_spread{i}; % LT spread
    idx_t                         = idx_t + add;
end

filename                            = sprintf('Matlab_panel_corrs.csv');
fullname                            = fullfile(PARMS.tabfolder, filename);
writematrix(panel, fullname)



